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Abstract 

Factorization of linear programming (LP) models enables a large portion of the LP 
tableau to be represented implicitly and generated from the remaining explicit part. Dy- 
namic factorization admits algebraic elements which change in dimension during the course 
of solution. A unifying mathematical framework for dynamic row factorization is presented 
with three algorithms which derive from different LP model row structures: generalized 
upper bound rows, pure network rows, and generalized network rows. Each of these struc- 
tures is a generalization of its predecessors, and each corresponding algorithm exhibits just 
enough additional richness to accommodate the structure at hand within the unified frame- 
work. Implementation and computational results are presented for a variety of real-world 
models. These results suggest that each of these algorithms is superior to the traditional, 
non-factorized approach, with the degree of improvement depending upon the size and qual- 
ity of the row factorization identified. 
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1 Introduction 



A recurring theme in the development of algorithms for linear programming has been the iden- 
tification and exploitation of special problem structure. Ideas as apparently disparate as the 
bounded-variable simplex method, primal and dual decomposition methods, pure and gener- 
alized network primal simplex algorithms, primal partitioning and column generation schemes 
may be unified to a degree with this view. 

The factorization approach introduced by Craves and McBride [1976] isolates special struc- 
ture in LP tableaus. We are interested in using factorization to reinterpret existing algorithms, 
and to discover common principles and apply them to develop new algorithms. Although all 
algorithms developed this way will, in theory, solve any LP, the efficiency of any particular fac- 
torization approach will be influenced by the relative number of factored constraints and their 
influence on the algorithm: the size and quality of the special structure isolated determines 
the influence of any particular factorization applied to any particular LP. 

Based on prior work by Brown and Graves [1975], in which generalized upper bound rows 
were successfully incorporated in a large-scale optimization system, we are interested in pursuing 
dynamic row factorization, where the dimension of the factored structure may vary (or even 
fail to be present) as the solution progresses. In our setting, we require the row structure of 
the model instance to be specified prior to solution, and that this structure remain fixed during 
solution. An extension of this approach is to allow the row structure to vary as the model is 
solved: this is a conceptually simple extension of the approach. 

Each algorithm is developed by factoring the constraints of the LP model into two classes: 
those that have the special structure (factored) and those that do not (explicit). This constraint 
factorization induces a factored structure in the LP tableaus which is exploited computationally. 
We demonstrate the dynamic factorization approach for three special structures: 

generalized upper bound rows; 

pure network rows; and 

generalized network rows. 

We implement each of the factorization algorithms by integrating it within the X-System 
(Brown and Craves [1975]). 

While the terms “partitioning” and “factorization” are frequently used interchangeably in 
the literature, we observe a distinction between the two approaches. We consider partitioning 
methods to be based on special structure in the original problem instance, which need not 
induce special structure in the LP tableau — in fact, the method need not be tableau- based. 
In contrast, factorization methods are based on special structure which occurs in bases and 
thus in the basic tableau. Thus, we classify dual decomposition (Dantzig and Wolfe [I960]), 
primal decomposition (Benders [1962]), and primal partitioning (Rosen [1964]) as examples of 
partitioning methods. 

Perhaps the earliest example of what we consider factorization is the treatment of simple 
upper bounds by Dantzig [1954] and [1963] and, independently, by Charnes and Lemke [1954]. 
They observe that it is more efficient to enforce the “logical” upper bound constraints with 
logical tests within the algorithm rather than treat them explicitly along with other “structural” 
constraints. While not originally presented in the context of a formal tableau factorization, the 
approach is easily viewed as such. 

The mutual primal-dual method of Graves [1965] focuses attention on the special role of 
nonnegativity constraints in linear programming. A clear distinction is drawn between the 
computational convenience of treating nonnegativity constraints implicitly rather than explic- 
itly and the unambiguous mathematical equivalence of all problem constraints, structural or 
nonnegativity. Emphasizing the special importance of inequality constraints, the approach 
yields an elegant theory and, as we will see, efficient implementations. We view this algorithm 
as the first formal example of factorization. 
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A similar primal-dual algorithm is presented by Balinski and Gomory [1965]. Related work, 
in which efforts are made to exclude slacks from the product-form representation of the primal 
basis, includes that of Zoutendijk [1970] and Powell [1975]. 

Dantzig and Van Slyke [1967] extend the earlier work for simple upper bounds and lend 
a more structured treatment to generalized upper bounds (GUB). In a problem with p GUB 
constraints and m structural constraints, their approach requires a working basis of dimension 
(m + 1), a considerable savings when p is large. 

Hartman and Lasdon [1972] specialize this approach to the multicommodity capacitated 
transshipment problem. In this case, the structure of the basic pure network columns in- 
troduces additional structure into the working basis, allowing further simplifications in basis 
representation and update techniques. Helgason and Kennington [1977] develop techniques for 
representing the working basis in product form and provide graphic interpretation of the basis 
updates. Kennington [1977] reports an implementation of the algorithm. 

McBride [1972] and Graves and McBride [1976] formalize and generalize the factorization 
approach. They view factorization as a unifying framework for tableau- based simplex special- 
izations and illustrate this by developing a variation of the GUB algorithm of Dantzig and Van 
Slyke and a GUB algorithm for the doubly-coupled linear programs of Hartman and Lasdon 
[1970]. They present a new algorithm for the set partitioning LP and an equality-constrained 
form of the pure network with side constraints model. Brown and Graves [1975] report an 
implementation of inequality-form, dynamic GUB row factorization for large-scale problems. 

Schrage [1975] extends the succession of simple and generalized upper bounds by introducing 
variable upper bounds (VUB), which are constraints of the form Xj < £*, where x * is said to be 
the variable upper bound of Xj . Schrage implicitly represents the VUB constraints by expressing 
VUB variables in terms of other variables. This permits the basis representation to be treated in 
two parts, one a large matrix which changes infrequently and thus needs only occasional update, 
and the other a small working basis which requires regular attention. Thus, computation and 
storage savings may be realized. Schrage [1978] extends these ideas to what he calls generalized 
VUB (GVUB) constraints, which arise frequently in models with fixed charges. 

Klingman and Russell [1975] sketch a factorization method for solving transportation prob- 
lems with side constraints. They suggest techniques for performing simplex iterations and 
updating the problem representation. Chen and Saigal [1977] present a similar approach for 
solving capacitated network flow problems with additional linear constraints. Both of these pre- 
sentations employ a graphical description of the basis update and treat the basis in two parts: 
one corresponding to a rooted spanning tree defined on the underlying graph, and the other 
a general working basis. Glover, et, al. [1978] report an implementation of the Klingman and 
Russell design, but one which (curiously) only accommodates a single side constraint. McBride 
[1989] reports an implementation which requires the pure network rows to be equalities and 
allows more than one side constraint. 

Generalized networks with side constraints are addressed by Hultz and Klingman [1976], 
who present details for the simplex priceout, column generation, and basis update. Hultz and 
Klingman [1978] report an implementation that (curiously) solves the “singularly constrained” 
generalized network problem. McBride [1989] reports an implementation that is not restricted 
to a single side constraint. 

The factorization approach has been extended by consideration of embedded structures. 
Glover and Klingman [1981] consider an LP with embedded pure network structure, i.e., the 
pure network structure appears in only a subset of the rows and columns of the technological 
coefficient matrix. They give an algorithm similar in spirit to their pure network with side 
constraint model, but the presence of the “side variables” significantly complicates the basis 
representation and update. They report an implementation of the algorithm but (curiously) 
restrict test problems to have no complicating variables. 

McBride [1985] solves an LP with embedded generalized network structure, presenting meth- 
ods for pricing, column generation, basis representation update and data structures. A success- 
ful implementation is reported to be about five times faster than MINOS (ca. 1977: Murtagh 
and Saunders [1977]) for the models tested. 
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Algorithms to solve problems with special substructures have motivated research to effi- 
ciently identify such substructures. Brearley, Mitra and Williams [1975] describe algorithms 
for detecting GUB row sets and exclusive-row structure sets (a set of rows whose structure 
may be transformed to GUB by column scaling). Greenberg and Rarick [1974] and Brown and 
Thomen [1980] develop algorithms to identify GUB sets. Brown and Wright [1984] identify 
pure network constraint substructures. Brown, McBride and Wood [1985] present a method for 
locating embedded and row-only generalized network structures. 

Todd [1983] develops a geometric interpretation of factorization which is for our purposes 
equivalent to the algebraic development of Graves and McBride [1976]. 

In the following sections we establish notational conventions, develop the mathematical 
foundation of a primal-dual simplex method and show the effects of row- factorization. Next, 
a general row-factorization algorithm is developed, specialized to GUB, pure network, and 
generalized network rows, and tasted on a suite of real-world problems. Finally, we discuss how 
the methods can be generalized further and how they can be applied more effectively. 



2 Mathematical Preliminaries 

The traditional statement of the linear programming (LP) problem is 

(LP) min : wy 

y 

s.t. a x y < r x , i = 1 , . . . ,m 

ejy>0 , j = 1, . ,.,n , 

where y is an n-vector of decision variables, w a vector of cost coefficients, each d{ an n-vector 
of technological transformation coefficients, each a scalar right-hand side coefficient, and ej 
the j th unit vector. While this statement of the problem is clear and unambiguous, there are 
reasons for preferring an alternative. The insistence upon drawing a formal distinction between 
the “structural” constraints aiy < r % and the “nonnegativity” constraints ejy > 0 obscures 
the mathematical structure of the problem by suggesting that the two types of constraints are 
inherently different. Certainly the exploitation of the special structure of the ejy > 0 constraints 
leads to computational efficiencies; however, in our theoretical development of the algorithm, 
we prefer to treat them simply as general inequality constraints. 

In order to achieve a consistent form, we rewrite the nonnegativity constraints as —ejy < 0 
and group them with the structural constraints. The problem statement then becomes 



(PLP) nun : wy 

s.t. : a x y < r x 



, i = 1, . . . , m + n, 



where wy is called the extremal function. 

From the standpoint of a primal algorithm, a matrix partitioned form of the primal tableau 
is derived. Let {aq, a X2 , . . . , a ln } be a basis for R n at y°. For notational convenience we will 
partition the constraints into two sets, those that are basic (binding) at y° and those that are 
nonbasic (not necessarily binding) at y° 
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Using this notation, the current basic solution y° may be expressed as By 0 = /, and since 
the rows of B are by definition linearly independent, y° exists. 

To isolate the important algebraic components let us assume that at the current basic 
solution the basis consists of h structural constraints and (n — h) nonnegativity constraints. 
Then 
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In partitioned matrix form, 
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and in partitioned matrix form 
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and we shall call DP the principal part of the tableau. By partitioning w = ( W\,W 2 ), g T = 
(<7iiP2) T and y° = (y°, y!,), the complete tableau may be written in partitioned matrix form as 
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Note that y\ is displayed explicitly in this tableau. Also, = 0 since the corresponding 
nonnegativity constraints are basic and thus binding. 

The corresponding dual problem is 



(DLP) max : xr 

X 

S.t. XQ? < tO* 7 , j = 1, . . . ,71 

xei <0 , i = 1 , . . . ,m . 

To develop a matrix partitioned form of the dual tableau, we proceed as before. Assuming 
the dual basis consists of h structural constraints and m — h nonnegativity constraints, we have 



T = = (a- 71 , a- 72 , . . . , a ? h , e ih+l , . . . , e im ) 

u = (u 1 , u 2 , . . . , u m ) = (to* 71 , to* 72 , . . . , to* 7 * 1 , 0, . . . , 0), 



so 



u = (u\ u 2 ) — (to 1 , 0). 



The nonbasic constraints are then 



K = (fc 1 , fc 2 , . . . , k n ) = (e il , e* 2 , . . . 

V = (n 1 ,!; 2 ,...,!; 71 ) = (0,0,... 



e %h , a- 7 ** 1 , . . . , a* 7 ” ) 
0, to* 7 ^ 1 , . . . , to- 7 ”), 
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SO V = (v l ,v 2 ) = (0, w 2 ) . 

The matrix partitioned form of the basis is then 
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and with the choice of Q = T 1 
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with the remaining constraints forming 
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The principal part of the dual tableau is 
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which we find to be exactly the principal part of the primal tableau, so 

DP = QK . 

Thus a single tableau representation supports both primal and dual algorithms. 



3 Interpretation of Primal and Dual Forms 

We may interpret a primal or dual algorithm as simply different perspectives of this same 
tableau, wherein a primal algorithm basis change is viewed as exchanging primal constraints and 
a dual basis change exchanges dual constraints. The classical Simplex Method may then 
be interpreted as solving (PLP) using the dual perspective. That the classical Simplex 
Method is naturally interpreted as a dual algorithm comes as a surprise to the conventionally 
trained. However, the consequent mathematical insight is compelling, especially in light of 
the notational simplification and apparent underlying role of A ^ , which we refer to as the 
transformation kernel 

There are several reasons for preferring a Primal-Dual algorithm to the Simplex Method. 
From a computational standpoint, because slack variables are carried logically rather than 
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introduced explicitly, we are able to clearly identify the essential information needed to execute 
the algorithm. The matrix A j”/ plays a key role in the calculation of the tableau, and the entire 
tableau can be constructed from Aj”/ and original problem data. Since Aj"/ is a submatrix 
of the inverse, T _1 , used by the Simplex Method, it is smaller and requires fewer arithmetic 
operations to update than does T” 1 . 

A second advantage of a Primal-Dual Algorithm lies in the flexibility it offers for special- 
ization to particular problem classes or structures. Indeed, it is the special structure and 
simplicity of the nonnegativity constraints that motivate the development of the algorithm in 
the first place. It is frequently the case that other special structures can be identified in classes 
of (PLP). Examples of such structures include simple upper bounds, generalized upper bounds, 
variable upper bounds, pure and generalized network substructures, etc. Such structure may be 
“static” in that its nature and dimension remains fixed throughout the solution process, or the 
structure may be “dynamic” in which case its precise nature and/or dimension may vary as the 
problem is solved. Some special structures may be more strongly characterized by their column 
structure and others by their row structure. The Primal-Dual perspective leads naturally to 
explanations of the implications of virtually any such problem structure and greatly simplifies 
the implementation of such a specialization. 

When an LP appears as a subproblem in a more sophisticated solution setting (for example, 
in a mixed integer programming problem or a nonlinear programming problem), the row/column 
symmetry of a Primal-Dual Algorithm is of critical importance in specializing the solution 
approach. The inherent symmetry of such an algorithm permits easy adaptation to branch- 
and-bound and cutting-plane approaches to mixed integer programming, to column generation 
settings, as well as to primal and dual decomposition techniques. 

We believe the reason for this flexibility offered by the algorithm lies in its more complete 
mathematical foundation. There is a natural consistency that arises from the choice of 
a vector space having the same dimension as the problem variables that is lacking 
in other approaches. A natural geometric interpretation of the solution trajectory follows di- 
rectly from this development. Incidental issues such as finding an initial basic feasible 
solution and dealing with degeneracy are resolved constructively in this math- 
ematical framework [Graves 1965], Other approaches resort to unnecessarily complicated 
tangential efforts. 

All the research results reported here can be developed, with some effort, in the framework 
of the classical Simplex Method. However, we choose to present these results in the manner of 
their development — the mutual Primal-Dual view presented by Graves. 



4 Column and Row Generation 



Rather than maintain a complete tableau DP, now consider the generation of just column c 
of this tableau. Rewriting in a manner that highlights our intentions, and labelling row and 
column partitions for identification 



DP = 
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\ — ^21 [^ 4 ] i 1 ] M2 — Aw] / 



By properly sequencing our computations we will exploit the fact that region ( it ) of a given 
column is simply a linear combination of terms in region (i) of the same column. 

Assume we want to place the current representation of column c into a work array z, which 
we partition as z T = (zf ,z%) to correspond to regions (i) and (ii). Expressed in terms of the 
transformation kernel AJ* , we compute column c as 
if c is in 
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and then 



21 = MhT . 



~2 = -j42i [ 21 ] ; 



or, if c is in ( jj ), 
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Then the current representation of column c is available in z T = 

The computation of row r of the tableau proceeds in a similar manner. We now view the 
principal part of the tableau as 
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If we want to place the current representation of row r in a work array z partitioned 
conformably with ( j ) and (jj) as z = (£3,24), we compute 
if row r is in (i), 



^3 — [^11 ]r > 

and then 

^4 = [-3]^ 12 ; 



or, if row r is in (ii), 



35 = [(-^Or^n'] , 



and 



z 4 ~ (^22)r + [-3]^12 , 



and the current representation of row r is available in z = (z.3,54). 

We see that in each case calculations proceed by first using a representation of to 
compute a portion of the row or column and then using this initial computation and original 
problem data to compute the remaining part. We will discover that our specializations ex- 
tend this approach by introducing additional tableau partitions which allow this computational 
strategy to be applied on a larger scale. 
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5 Transformation Kernel Update 

The dynamic behavior of A^ 1 is important. We see from the primal row basis B and non basic 
rows D that the dimension of A^ 1 corresponds to the number of basic structural constraints, 
or, equivalently, to the number of nonbasic nonnegativity constraints (recall that if a nonnega- 
tivity constraint is nonbasic and thus nonbinding, the corresponding variable may possibly be 
nonzero). Recalling that our primal view of a basis exchange is as an exchange of constraints 
between B and D, we see that one of four cases may occur during a pivot 

A structural constraint enters the basis B and a structural constraint leaves the basis and 
enters D. Since the number of basic structural constraints (and the number of nonbasic 
nonnegativity constraints) remains unchanged, the dimension of A is unchanged. A 
pivot of this type involves a row in region (j) of B and a row in region ( ii ) of D, and thus 
it corresponds to a pivot coordinate in the location ((ii), (j)) of the tableau DP. 

A nonnegativity constraint enters the basis and a nonnegativity constraint leaves the 
basis. Again, the dimension of A J*/ remains unchanged. Since this pivot involves a row 
in region (jj) of B and a row in region (i) of D, the corresponding tableau DP pivot 
coordinate lies in ((i), {jj))* 

A structural constraint enters the basis and a nonnegativity constraint leaves the basis, 
and thus the number of basic structural constraints (equivalently, the number of nonbasic 
nonnegativity constraints) increases by one . The dimension of A ^ is increased by one. 
This corresponds to a pivot coordinate in region ((H), (jj)) of the tableau DP. 

A nonnegativity constraint enters the basis and a structural constraint leaves the basis, 
and thus the dimension of A J"/ decreases by one . The corresponding pivot coordinate in 

DP is ((i),0)). 

We see that we may exert some influence on the behavior of the dimension of A^ 1 by our 
strategy for selecting target exchanges for primal and dual constraints (i.e., our pricing strategy) 
and through our tie-breaking rules for choosing pivot row/column, and that this dynamism is 
an inherent feature of an effective algorithm. We have already seen the fundamental importance 
of the kernel (An) in our computations. Thus, a successful implementation must manage this 
dynamic behavior efficiently and reliably. 



6 Factorization 

The row- factorized problem to be considered is 
(FLP) min : wy 

s.t. : Ey < r } explicit constraints 

Fy < b } factored constraints 
-ly < 0 } nonnegativity constraints , 

where y is an n-vector of decision variables, w a vector of cost coefficients, E a matrix of con- 
straint coefficients for “explicit” constraints with right-hand side m-vector r, F a matrix of 
constraint coefficients for “factored” constraints with right-hand side p - vector 6, and — / the 
negative of the identity matrix. In this general development, we refer to the F-type constraints 
as “factored” only to distinguish them from the “explicit” F-type constraints, and assume noth- 
ing about their structure. Not until our specializations later will we impose special structure on 
F, and the structures we will consider may permit the representation of the F-type constraints 
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without the inversion of a matrix. We will see that this approach is centered around handling 
the part of the basis corresponding to the £-type constraints explicitly while factoring the por- 
tion of the basis corresponding to the T-type constraints. The notation is chosen to suggest 
this idea. 

Recall that a basis for the primal algorithm consists of n linearly independent rows from 
the constraint matrix when it is assumed to include both structural (explicit and factored) and 
nonnegativity constraints. Assume that the current row basis consists of k rows from £, l rows 
from F and (n — (k + /)) rows from — /. Repeating our notation 



B = 



k+l 

f 

An 

0 



n—(k+l) 

\ 

A \2 } k + 1 

-/ J }n — (k + l) 



where [A\\ A12] includes all basic structural rows, from both E and F. 

We will ultimately be interested in isolating the effect of each type of structural constraint 
algebraically in the factored tableau, and thus we require greater resolution in our factored 
basis. Introducing obvious notation, we have 



/ 

[An A\o] = 

where the kernel of dimension ( k + l ) is given by 



n-(k+l) 



E u 


E 12 


E 13 




}k 


Fo\ 


E 22 


E 23 


) 





[^11] 



E 11 £12 

F ‘2 ] F22 



Because An is a basis for R k+l it follows that it is always possible to identify among the 
columns of [F21 F22] a nonsingular submatrix F02 of dimension /, since otherwise the rank of 
[F21 ^22] is at most (/ — 1) and thus the rank of A\\ is at most (k + l — 1), or equivalently 
the rank of B is at most (n — 1). We will later see that one of the important implementation 
challenges is the task of efficiently managing the structure and nonsingularity of F22* 

The full factored row basis is then 





k 


1 


n-(k+l) 




(" 






U) 


E 11 


E 12 


E 13 


(jj) 


T21 


E 22 


E 23 



( 333 ) \ 0 0 



-/ 




(k + /) 



Introducing the notation 
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An = En — E\2F 22 l £21 
Al3 = £l3 — Ei2F 22 F 23 , 

where yin is the Schur complement, or Clauss transform, of F22 in A\\ (e.g., Golub and Van 
Loan [ 1983 ]), we can write its inverse as 



B ~ 1 



A 1 / — E\2F 22 A \3 

-FgFnAj (I + F 22 'FnA^En)F 22 ] F 22 \F2z - FuAjAn) 

0 0 -I 



Grouping the coefficients of the nonbasic constraints and applying the same column ordering 
yields 



k l n-(k+l) 



(0 


-I 


0 


0 


}k 


(™) 


0 


-I 


0 


}l 


(Hi) 


£31 


£32 


£33 


} m — k 


{iv) 


y £41 


£ 42 


£43 J 


}p-i 



The principal part of the factored tableau is DP , where P = —B 1 is the conjugate row 
basis. With the additional notation 

A3] = £31 - Es2F 22 F21 

A33 — £33 — £32 F^ 1 F23 

Ai = £41 — F42F22 ^21 
£43 = £43 - Fi2F 22 F23 , 

the principal part of the factored tableau is 





O') 


0^ 


Oii) 




( 




\ 


W 


4-1 

yi n 


-A^E] 2 F 22 l 


in 1 ^13 


DP = 


-£ 2 - 2 , £ 2 iiri 1 


(/ + £ 2 - 2 1 £ 21 iri 1 £i2)£ 2 "2 1 


£22 ' (£23 — £21 -An A13) 


(Hi) 


-j 43 ij 4 n 


(A3i J 4 ] f 1 1 £i2 — Ez 2 )F 22 x 


^33 - is.ir.'iis 


(iv) 


^ -Ai^n 


(£41 ^ 4 ] i’ £12 - £ 4 2 )£ 22 1 


£43 — £41^7/ j 4 i 3 J 
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Partitioning 


W = (nq, w 2 ,w^w 4 ), r T 


— (rf,r 2 ) and (F — 


(bf,b 2 ) and introducing 


notation 










u>2 = 


w 2 - w\F 22 F 21 






U)z = 


- w 2 F 22 F 2 z 






62 = 


b 2 — F41 F 22 b\ 






f\ = 


i' 1 — FiiF^&i 






*2 = 


r 2 — F 2 iF 22 1 6] , 




the complete factored tableau is 






r 

n \\ 


-A u ' E\ 2 F 22 




Mi'n 1 


—F22F21-A.il 


(/ + F 2 - 2 l F 21 i 1 - 1 1 E l2 )F 22 1 


C 2 2 1 _(F 23 — F 2 \A^[A\z 


1) f 22 \f -f 2 } a u 'f) 


-A 3 \Ail 


(^ 3 Mh'F] 2 - F 32 )F 2 - 2 i 


Mi — M\A~^i A\z 


r 2 - As\Ai\V\ 


-Fa Mu' 


( F41 F 12 — F 42 )F 22 ' 


C43 - F 4 )j 4 j'j 1 ^ 4 i 3 


~h - F^Aii'i-] 


-w 2 Aii 


(w 2 AuE 12 - W] )F 22 ' 


wz - w 2 A^A n 


w\F 22 b\ + w 2 Ai |’fi . 



We see that with knowledge of the current factorization, we can construct the entire tableau 
from F22 > anc * the original problem data. The dimension of F 2 2 is equal to the number of 
F-type constraints that are currently basic, and thus can be at most p. The dimension of A u ] 
is equal to the number of F-type constraints that are currently basic, and thus cannot exceed 
m. We call Af' the explicit transformation kernel and f 22 the factored transformation 
kernel 



7 Factored Column and Row Generation 

Consider generation of column c from the principal part of the tableau DP . Rewriting in a 
manner that highlights our intentions 





0 ) 


Ui) 


Un) 


(0 




[-A^EnF^} 


_ - \ 

[^11 1j4 13] 


(«) 


{-f 22 ] f 21 [ ]} 


{f 22 — F 22 F 2 \[ ]} 


{F 22 'F 23 - F 2 - 2 ‘F 21 [ ]} 


(m) 


_ ^31 [ ] — F 32 { } 


— F 3 ][ ] — E 22 { } 


F 33 — F 31 [ ] — f 32 { } 


(in) 


-F 41 [ ]-f 42 { } 


-c 4 i[ ]-F 42 { } 


c 43 — F 4 i[ ) - F 42 { } J 



where A \ 3 is defined as before, and the brackets “[ ]” and “{ }” contain terms common to but 
displayed only once for each column. 

Assume we want to place column c into work array z. We partition z conformably as 
z T = (zj , Zo , zj, zj), refer similarly to components of unit vector e c , and employ a (m + p)- 
vector work array z. The notation ” denotes simple assignment, indicates that a set of 
factored equations must be solved, and “=>” explains the corresponding result. 

If column c is in (j), 
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and then 



- [ir.t, 



then 



z 4 F2\[z\] 

F22Z2 = z => Z2 { — ^ 7 22 1 ^ r 21 [-1] 



Z3 4 ^31 [z\] - Es 2 {z 2 }, 



and finally 



if c is in (jj) y solve 



24 < — £41(21] 


- £42(22}; 


F 22 z 2 = -(e 2 ) c 




22 - -(£ 22 1 ) c 


z E\ 2 z 2 


=> 


z < (E 12 £’ 22 i ) c 


z\ «- A^z 




21 - [-iri 1 ^2£]T2 , ] c 



then 



and finally 



2 <- (e 2 ) c - F 2 i [zi] 

F 22 z 2 = z =» z 2 - {(F 2 - 2 ') c - F22 F 2 \ [zj]} , 



Z3 < £31(21] — £32(22} 

24 * £41(21] - £42(22}; 



if c is in (jjj), 



2 - (F 23 ) C 






£22 2 2 = Z 




22 - F 22 1 (F 23 ) C 


2 <— (£l 3 ) C ~ E\ 2 Z 2 


= 4 > 


2 ( z'l 13 ) C 


21 <— ^n '2 




21 <- [^ndil 3 ] C , 



then 



2 <— (£ 23 ) c - £2121 
F 22 z 2 = z 



^2 ^ {F 22 \F 2Z y - F^F 2 ,[z x }}, 



and finally 
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-3 <— (F 33 ) c - Ez\[z\] - E32 { ~ 2 } 

^ 4 -(F 4 3 ) C -F 4 ,[ri]^£; 42 {^}. 

Similarly, row r can be generated. The principal part of the tableau is now viewed as 





U) 


(ii) 


(iii) 


(i) 


[^n'l 


{-[ 


\ 

[ JF .3 + { } F23 


{ii) 


l-Fv'FnA^} 


| {F 2 - 2 '-[ JF.2F22 1 } 


[ JF .3 + { }F 23 


(iii) 


[—431 ^n 1 ] 


{(-F 32 -[ ]F 12 )F 2 2 1 } 


F33 + [ ]F ,3 + { }F 23 


(iv) 


\ [-AMfi 1 ] 


{(-F42 - [ ]F,2)F 22 1 } 


F43 + [ JF.3 + { } F23 J 



where A 3 i and F 4 i are defined as before, and the brackets “[ ]” and “{ }” contain terms common 
to but displayed only once for each row. 

Assume we want to place row r into work array z. We partition £ conformably as z — 
(£5, i?6, z?), refer similarly to components of unit vector e r , and employ a (m + p)-vector work 
array z. 

If row r is in (i), 



K 1 1 ] 



then 



z * z§E\2 =$> z * [A^rEn 

Z&F22 =Z => Zfi ‘ {-[A u l }rE u F 22 '}, 



and finally 



h *— [h\E \2 + {SbJF*; 



if r is in ( ii ), solve 



then 



and finally 



Z6F22 = -( e e)r 
Z <— ^ 6-^21 

~5 — zA u l 



Ze <- -(F£)r 
Z - -(F^) r F 2l 
Z5 [-F22* F2 l^ll'lr, 



z ■*— (ee) r - [zs\E\2 

Z6F22 = Z => Z 6 *- { (^*22' )r - [z5)E n Fn ] }, 
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[-= 5)^13 + { -^6 } ^ 23 ; 



Z 7 <- 



if 7- is in (Hi), 



then 



Z * (-^32 )r 

h F 22 = z 

Z <— (£ 31 ) 1 - + -?6^21 

55 - Mr, 1 ] 



ie - -(£32)^22' 

Z *— (A 3 I )r 

h [-^ 31 -^n]r 



5 - -(£?32)r - [5 s]£?12 

Z&F 22 = z => Z 6 <— {( (-£'32)r - } 



and finally 



z-i <— (^33) r + [^5] -Ei 3 + { 5 6 }F 2 3 ; 



if 7- is in (in), 



then 



Z <■ (i ? 42)r 

Z6^22 = Z 
Z <— (Fn) r — Z6 F 2 i 

i 5 <- \zA^\ 



=> h <— -(F 4 2) r F 22 1 

=> ^ (Al)r 

£5 * [ -^41 ^1 1 jr 



Z * (^42 )r — [zs\E\2 

^ 6^22 = Z => Z 6 <— {(-(F 42 ) r — [^5 ] ^1 2 ) ^22 1 } ’ 



and finally 



£7 (F43) r + [z$]E \2 + {ze}F 23 - 



8 The Complete Algorithm 

The complete algorithm is described in terms of abstract functions which operate on funda- 
mental data structures. 

Tableau management requires two index maps: one yielding the intrinsic coordinate in the 
principal part of the tableau for each original, extrinsic problem row or column, and the other 
its inverse map. Intrinsic arguments are shown in lower case, and extrinsic in upper case. 
Index_Exchange(indexl,index2) updates these maps for the exchange of a pair of tableau 
coordinates index 1 and index2. 

The tableau regions are successive partitions of indices 
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TABLEAU 

REGION 


ENDING 

INDEX 


CONTENTS OF REGION 


(0 


MEC 


basic Columns solving Explicit rows 


(ii) 


MFC 


basic Columns solving Factored rows 


(Hi) 


MER 


nonbasic Explicit Rows 


(iv) 


Tit 


nonbasic Factored Rows 


U) 


NER 


basic Explicit Rows 


(ii) 


NFR 


basic Factored Rows 


(Hi) 


m + n 


nonbasic Columns. 



Increment(endingJndex) and Decrement(endingJndex) are functions to modify these 
ending indices. 

Generate_Row(row) and Generate_Column(column) place numeric values of a tableau 
row or column in ROWCOLQ, which is commensurate with the tableau dimensions. 

Using ROWCOL(), Update_Rim(row,col) maintains current numeric values of the right- 
hand-side and bottom row of the tableau in RIM(). 

The explicit transformation kernel, is operated on by functions using ROWCOLQ: 

Add_E-Kernel JRow(ROW), 

Add_E-Kernel_Column(COLUMN), 

Delete_E-Kernel_Row(ROW), 

Delete_E-Kernel_Column(COLUMN), 

ReplaceJE-KerneLRow(REPLACEDJtOW, REPLACING JlOW), and 
Update JExplicitJIVansformation_Kernel, the pivotal update. 

A ^ can be represented in any way that suits the implementer and efficiently supports these 
functions. Generally, we find to be relatively sparse, and report here results using an array 
with an entry-point for each inverse row and column which accesses a stack of orthogonally- 
linked nonzero inverse elements. We have also implemented a dense vector-processor version, 
and the design issues for LU-based schemes are given by Olson [1989]. 

The factored kernel, F 22 , is manipulated with: 

Factored_KerneLSingular(ROW, COLUMN), a Boolean function, 
Find_E-KerneLColumn_for_Key (ROW, COLUMN), a column from region (i), 
Find-F-Kernel_Column_to_Remove(ROW, COLUMN), a column from region (ii), and 
Update_Factored_Kernel(ROW, COLUMN), the pivotal update. 
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The complete abstract algorithm is: 



0. Initialize; 

1. select primal or dual algorithm; 

2. Select a primal (col) or dual (row) violation; 

STOP if current solution is terminal, or 
Generate_Column(col) or Generate_Row(row); 

3. By ratio test with RIM() and ROWCOL(), 
select a (row) or (col) pivot coordinate; 

STOP if current solution is terminal, or 
Generate_Row(row) or Generate_Column(col); 

4. Update_Rim(row,col), 

primary Index_Exchange(row,col), 

perform secondary and tertiary exchanges (Table 8-1), 

Update_Factored_Kernel(ROW,COL) 

Update JExplicit_TVansformation_Kernel, 

Go to step 1. 



18 



Oii) 



S' o 

<3 o 



s 

o 

K 

j 

o 

o 



D 

O 

sc 

* 

o 

* j 

5 I 
® w 
a ® 



o, ^ o q 

asaf 

JS3|*j 

w "a ® u 

x E itf ® T x 

4 Si g sh 

x £ • ® « x 
® u *® 3 * 3 . ® 

"3 Jf -all's 

— a u, o cc £ 



- .a 

£9 



-S, 



S j££ 



® s a c 2 2 



11 



n 

1 1 

l if! 

aa 



i 

G 

-S 



51 

g§„ 

3= fe | 

g|S 

w n z 

?iE~ 

iiii 
*1453. 
13 S a 
ssis 

* SSI 

*2 * h x 

® “ * X . 

“ •* *3. ® ^ 

fifilSl 

fc u 



a o 

=a 

o 



o 

§51 

J o 

o 2 

o ® ^ «*• 

7 £ 



S' 

o 

sc 

z" 



*44 



o: 

Lb 

z 



^ia 

» at g ■» 
icu. ® o 

2 5.^ S 
5s2¥ 



3 I ||. gg 

’ 4 ? 5 I 

n J « 

H w ** * 



I a £ . 

jih 

j t£S 

ill" 



« > I q 
» ®i ui “8 



^.o 

*2 

s, s,£ £ 

Jlrr 

« » c a 



4 2 8^ x x = = 

3 sISu-h^ 1 3 

•9.11*35 5 8 e e 
u 0 ae ° £ -S -S g £ 

«QQ 



o 

p 



C* 

£zg 

z cd S 



1 _ -2-5- 

3s£ ||5'‘" 

S55 g 2 ! 

“ — j 2 

® tq *G w U 9 

c a M M V 

* ® ® cj C*1 * 

4 §§ 11 ^ 

ii ill 1 



o' g 

y 2 
^.o 

.jt 

Up * ® 
U H 5) » 

3 2 2 3 
r~gg 

c a n h 

1111 
® ® g g 

tj C -o *u 
a a a a 



s 

o 

* 

j 

o 

u 



0 ~ a 2 

a sr4 

1 “J 

3 "3 j 3 
o ^ o 

u£.4 

Hi! 

® <** U ® 

¥ g cS ^ 

s g|I 



?p 



2 

■? 

|I 

56? 

8 IS 

“4 * 

.. siE~ 
s i] * i 

mu 

^IHe- 

£ | J £ * ! 

z. s* » ti ~ 
“J'S* . 5 
d Si Sii 84 

• siihE 

jS2ff 5 



r 



aT ul 
» u 



p z « ® 

1 CJ U M 

S5.S S 



Sigg 

s ! « i 

U U *U *U 



*T i ® o 

fjtgii 

■? !l?St 

• x « S « S 

® ^ 2 i it 



33$ 

H 5 

?o 

oO 

5 ? 

j o 
O a 

u S 
X°i 

4 ^ 

3 a 

» a 

S 3 
*1 

13 

® e 
^ § 

oh 

L> -O 

AS 



o 

>»o 

AT- 



^■S ? 
9*o 

JSi 

3 o -i O 
i n c 

C ® 

V — k t£ 

» « ® a 

i o c ^ e 

; «i sS-e 

' ® SC * x 

• si 8H 

® i J X . 
a "O *5 e u. 
»2 «*Oq 

u 



(A 

<D 

GO 

C 

cC 

J2 

u 

X 

W 

d 

a 

JH 

3 



& 

'O 

c 

CO 

Sm 

CO 

•n 

c 

o 

u 

<D 

tfl 



^ oj ^ 

CO — Q) 

u ^ u 

ai a o 



3 < 



o 

£ 



!g| 



i as 

it™ 

* ® z £X 

M HU Lb 



o£ 

y s 
o 
» u 



5 M MU Lb 2 

2 2 2 5, z 2 

^ ^ J3 ^ 

5 x 8 § 

“ a 8“ 

Js X X « 9 Z. 
2. ® ® 5 U J> 
® ^ ^ ® ® ® 
Q»QQQ 



S 3 

J= J3 
u u 



2 2 



8 8 
£ £ 
8 8 
a a 



0 

£^2“ 
9<ty 
9 u ~ 

3 z a: 

uli,_ 

1 &&£ 

§ s 9S 

^ J3 J3 ^ 
* W U B 

ii x x ; 

111 
® ® ® r 

®‘2'2 ^ 
q ° * a 



z 

2 

§ 

a: 



Ills 

£S g S 

O c 

!««!! 

^ II || ^ Sd 

QkQiiCQ 



8 

2? 




ao 

W 

n 

ss 



19 



aod Update-Explicit-Trausformation-Kernel. 



Functions for maintaining the factored kernel vary with the factorization, and a good im- 
plementation will exploit these differences. However, a general specification will suffice for all 
factorizations here. 

We are exclusively interested in performing two fundamental operations: 

1. Solving factored equations of the form: 

F22-2 = k >2 

and 

%F, 22 = 62 

where Z2 and Z2 are unknown and 62 and 62 are rational (not necessarily integer), and 

2. Restoring F22 to the desired form of F22 which makes the factored equations above easy 
to solve, where F22 results from inflicting F22 with a column exchange, a column and row 
deletion, a column and row addition, or a row exchange. 

Factored _KerneLSingular(LABELl,LABEL2) predicts whether F22 will be singular if: 

1. LABEL1 gives a basic column in a basic factored ROW for which COL=LABEL 2 would 
be exchanged; 

2 . LABEL1 is a basic factored ROW solved by basic column COL=LABEL 2 , both of which 
would be removed from the factored kernel; 

3 . LABEL 1 is a nonbasic factored ROW, which would be added to the factored kernel with 
column COL=LABEL2, or 

4 . LABEL1 is a basic factored ROW which would be replaced by some other row LABEL2. 

The first case, a column exchange, is equivalent to asking whether solving F22Z2 = 62 with 62 
equal to region ( ii ) of the proposed entering column COL, yields a nonzero term associated with 
the basic factored row in ^(row). This is easy to answer for the factorizations we discuss. For 
instance, Brown and McBride [ 1984 ] show that back-solving a “nearly-triangulated generalized 
network” basis F22 only as far as ROW will suffice, and that this can be implemented as a 
traversal of the “backpaths” of the (zero, one, or two) coefficients in b 2 for the column now 
basic in ROW. If numerical precision is an issue, ^(row) can be computed during this search 
and tested for significance. 

The second case, a row and column deletion, is trivial because the resulting F22 will always 
be nonsingular. 

The third case, a row and column addition, can be answered by adding ROW and COL to 
F22 without disturbing its desirable special structure, thus creating an instance of case one. For 
instance, placing ROW first, and COL last in F22 suffices for all factorizations discussed here. 

The fourth case, a row exchange, can be answered by applying the second case, deleting 
ROW and its basic column, and then (perhaps repeatedly) applying the third case, adding the 
row with index COL and any column which would yield a nonsingular F22. 

Find-F-KerneLColumn_to _Remove(ROW,COL) given basic factored row ROW, re- 
turns its basic, or “key” column COL. 

Find_E-KerneLColumn Jor_Key(LABEL,COL), given either a nonbasic factored ROW= 
LABEL, or a basic column LABEL solving a basic factored ROW, searches for an acceptable ba- 
sic column COL in Explicit Rows (region (i)) using Factored_Kernel_Singular(ROW,COL). 

Update_Factored_Kernel(LABELl,LABEL2) restores F \ 22 to the desired form of F22, 
with a possible increase or decrease in dimension. Following the case- by-case scheme of Fac- 
tored_KerneLSingular, we pre- and post-process the factored basis representation to permit 
use of a single, static factorization update of conventional design. Olson [ 1989 ] pursues this in 
considerable detail. 

Table 8-2 displays supporting data structures for these functions. 
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TABLE 8-2 Factorization Algorithm Data Structures 



FACTORIZATIONS 


DATA 

STRUCTURE 


SIZE 


USE 


GUB,PN,GN 


RIM() 


m + n 


current tableau right-hand side, 




ROWCOL() 


m + n 


bottom row 

current tableau row and column 




MSKRC() 


m + n 


logical mask true for corrasponding 




LQRC() 


m + n 


nonzero in ROWCOL 0 , false otherwise 
LIFO queues of nonzero row and 




KEY() 


m + n 


column coordinates in ROWCOL() 
basic column in basic factored row, 


PN,GN 


WORK() 


V 


and vice versa 
values for 62 or bo 




MSKWK() 


P 


in basic factored equations 
logical mask for nonzero row and 




LQWK() 


V 


column coordinates in WORK() 
LIFO queue of nonzero coordinates 




PO() 


V 


in WORK() 

next basic factored row in pre-order 




P() 


V 


off-diagonal row with 




D() 


V 


nonzero factored coefficient 
depth, remaining back-substitution 


GN 


VGN() 


p 


path length in factored component 
generalized network cycle factors 




JMUL() 


n 


ratio of generalized network 



coefficients in each column 



Generalized upper bound (GUB), pure network (PN), and generalized network 
(GN) factorizations respectively require more data structures to support kernel 
factorization. For each factorization, F 2 2 is maintained in some partial ordering 
of rows, and of columns — a signed identity for GUB, upper-triangular for PN 
(e.g., Bradley, Brown, and Graves |1977j), and nearly-upper-triangular for GN (e.g., 
Brown and McBride [1984]). Direct solution of factored kernel equations F77Z7 — b 2 
and zj F 77 = is performed alternately using the data structures shown. 
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9 Computational Experience 

The factorization methods introduced here have been implemented and used to solve a variety 
of existing models provided by our colleagues. With their help, we have extracted suggested 
problem instances from a diversity of decision support systems on host computers ranging 
from mainframes to micro-computers. Because our goal is to test factorization technology in 
isolation, the results reported here are achieved without benefit of any model-specific knowledge 
or tuning. 

However, we also seek to develop effective modeling tools for customized use in developing 
and refining new models. To this end, we have greatly benefited from the experience and advice 
of our colleagues, and we include some discussion of modeler guidance and insight along with 
the numerical results. 

Each model is introduced below by a short synopsis. Multiple instances of some models 
are reported where diversity of size, structure, and taxonomy have proven interesting to the 
modelers. For those models that employ nonlinear, mixed integer, or decomposition features, 
we report solution statistics for the initial linear program. 

• GTE The seven Telephone Operating Companies within GTE have adopted an integrated business system 
called Capital Program Management System (CPMS) to guide their 3 billion dollar per year capital planning. 
The system includes a large scale mixed integer programming optimization system that optimizes the critical 
economic tradeoffs between maximizing the long-term budget value of the firm’s equity and satisfying shorter- 
term financial constraints, resource limitations and service objectives. Investment opportunities for the next 
5 years are modeled as 0-1 variables with alternative implementations for each. The objective is to maximize 
the net present value of the capital portfolio. There are financial constraints on capital, internally generated 
funds, net income to common, and limits on resources such as labor hours, lines installed, etc. There are also 
constraints that enforce logical relationships among opportunities (such as, if choose A then must choose B). 
See Bradley [1986]. 

• INVEST Capital allocation and project selection for Mobil Oil Corporation are modeled as a two-stage 
multi-year nonlinear capital budgeting problem with over 40,000 integer variables. A master problem alio* 
cates capital among markets over a multi-year horizon considering the estimated nonlinear effects on sales 
of concentrated marketing investments. The instance reported here is a mixed- integer linear program sub- 
problem of the two-stage model which, given these annual capital expenditure limits for a market, selects 
particular alternate investments. Such subproblems are easy to solve, and optimality is achieved with a single 
iteration of the nonlinear master problem. See Harrison, Bradley and Brown [1989]. 

• TANKER A crude oil tanker scheduling problem faced by a major oil company is solved using an elastic 
set partitioning model. The model takes into account all fleet cost components, including the opportunity 
cost of ship time, port and canal charges, demurrage and bunker fuel. The model determines optimal speeds 
for the ships and the best routing of ballast (empty) legs, as well as which cargoes to load on controlled ships 
and which to spot charter. All feasible schedules are generated, the cost of each is accurately determined 
and the best set of schedules is selected. See Brown, Graves and Ronen [1987]. 

• HFDF A large-scale elastic set partitioning model used to assign frequencies for a network of high frequency 
direction finding receivers. See Brown, Drake, Marsh, and Washburn [1990]. 

• GAS A multi-time period strategic model for use by natural gas utilities for determining optimal contract 
levels for gas purchase, storage and transmission. An underlying generalized network flow model represents 
gas being bought, stored, shipped and consumed over a multi-year time horizon, typically at a monthly level 
of detail. Constraints and variables are added to handle variable maximum and minimum purchase levels, 
variable leased or constructed storage and variable transmission capacities. An integrated parallel model 
incorporates the peak requirements necessary on some days during cold winter months. This model has been 
used by a number of utilities including Southwest Gas Corporation and Questar Pipeline Corporation to plan 
operations and to justify such plans to regulatory agencies. See Avery, Brown, Rosenkranz and Wood [1992]. 

• KELLOGG A multi-time period, multi-plant production/inventory/transshipment linear program for Kellogg 
cereals. The model guides weekly processing, packaging and shipping decisions. Production consists of two 
stages: processing lines produce basic products which are then packaged on packaging lines into different- 
sized containers to yield finished products. Processing lines produce a subset of the basic products and have 
limited capacity with overtime charges for weekend shifts. Packaging lines are analogous. In-house inventory 
capacity is limited although outside storage is available at additional cost. Inter-plant shipments of finished 
products are made by rail or truck. See Wood [1989]. 



22 



• ODS A commonly occurring problem in distribution system design is the optimal location of intermediate 
distribution facilities between plants and customers. A multicommodity capacitated single-period version 
of this problem is formulated as a mixed integer linear program. A solution technique based on Benders 
Decomposition is developed. . . . An essentially optimal solution is found and proven with a surprisingly small 
number of Benders cuts. See Geoffrion and Graves [1974]. The instances reported here are decomposition 
master problems. 

• TAM The annual decision on how much the Air Force should spend on aircraft and on munitions is of great 
interest to many people. How the Air Force staff develops information to support the decision has changed 
over the years. Currently, a linear program is being used by the Air Force Center for Studies and Analysis 
and is being tested by the Munitions Division of the Plans and Operations Directorate (AF/XOXFM) for 
munitions tradeoff analysis. The LP uses existing data and estimates of (1) aircraft and munition effectiveness, 
(2) target value, (3) attrition, (4) aircraft and munition costs, and (5) existing inventories of aircraft and 
munitions. Other factors considered are weather and length of the conflict. See Might [1987] and Jackson 
[1989]. 

• PHOENIX A planning model for the multi-year, multi-billion dollar modernization of the U. S. Army’s aging 
helicopter fleet. The mixed integer linear program employs a multi-product production/inventory formulation 
with aged inventory. Goal constraints attempt to enforce fleet size, maximum age, and technology goals for 
each year and each of four aircraft missions, while also keeping expenditures within upper and lower limits. 
Additionally, combinatorial constraints and variables handle production line startup and shutdown costs, 
minimum and maximum production levels and requirements linking certain production lines. See Clemence, 
Teufert, Brown and Wood [1988] and Brown, Clemence, Teufert, and Wood [1990]. 

• EA6B Configures jammers of hostile radar on an EA-6B “Prowler” Naval electronic warfare aircraft. See 
Sterling [1990]. 

• DEC Digital Equipment Corporation uses this model to determine worldwide manufacturing and distribution 
strategy for new products. This mixed integer, linear program suggests a production, distribution, and vendor 
network which minimizes cost and/or cumulative cycle times subject to constraints on estimated demand, 
local content, and joint capacity, over multiple products, echelons, and time periods. Cost factors include 
fixed and variable production charges, distribution via multiple modes, taxes, duties and duty drawback, and 
inventory charges. See Harrison, Arntzen, and Brown [1992]. 

• AMMO 4H A four-commodity transshipment model for delivery over time of military products from pro- 
duction and storage locations to overseas locations to support theater operations is developed. The model 
covers five physical echelons, including production plants, storage depots, ports of embarkation, ports of 
debarkation and geographic field locations. Road, rail, sea and air transportation are modeled, and product 
demands are time-phased. Capacitation occurs primarily on sea and air links, and as throughput capacities 
on transfer points, requiring replication of some echelons. The objective of the model is to minimize deviation 
from on-time deliveries. See Staniec [1984]. 

• BUSCH A model of brewery- to- wholesaler movements of beer for Anheuser-Busch. The model also includes 
some packaging decisions and is essentially a multicommodity flow model with joint capacity constraints 
arising from loading dock and inventory capacities as well as some managerial requirements. See Brown, 
Mamer, McBride and Wood [1992]. The instance reported here is a small pilot model for the full-scale 
system with millions of variables which is solved directly, or by decomposition. 

• BAR A linear, mixed-integer multi-period production-inventory master planning model. See Harrison [1992]. 

Four implementations are compared: “XS” is an unadorned version of the X-System, an 
implementation of the Graves mutual primal-dual method with its GUB factorization disabled, 
while “XS(GUB)”, “XS(PN)” , and “XS(GN)” each employ the respective factorizations dis- 
cussed here. To establish a frame of reference, performance of these implementations is com- 
pared with two well-known commercial solvers: IBM’s Optimization Subroutine Library “OSL” 
(Release 2 [1991]), and two versions of “CPLEX” (Version 1.2 [1990] and Version 2.0 [1992]). 

Ideally, one would develop four equivalent formulations of each model, each customized for 
its particular solver with the goal of inducing a large factored row set of the appropriate type. 
This approach is a consistent theme in the literature dealing with specialized algorithms and 
one that we strongly endorse. Alternate formulations of a model are often available, and it 
seems sensible to choose one that exploits as much as possible the strengths of the solver. 
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However, all of the models used here are “off-the-shelf” in the sense that they were developed 
at various times by various modelers, and alternate formulations are impracticable. Thus, the 
approach is to preserve a single, unfactored representation of each model, and attempt to 
identify favorable row structures through the use of heuristics. The procedure is based on 
the work of Brown and Thomen [1980], Brown and Wright [1983] and Brown, McBride and 
Wood [1985]. The heuristics are greedy and myopic in the sense that they initially consider the 
entire row set of the problem, and discard one row at a time without backtracking until the 
remaining set satisfies the desired row factorization. This can be expected to confound or destroy 
structure introduced by the modeler. Although the automatic factorization implementation has 
options to accept modeler guidance, the methods are compared here without this subjective 
complication. While these model-naive experiments yield interesting and useful observations 
about the implementations, they suffer for lack of guidance by a skilled modeler. 



24 



Table 9-1 shows the important structural information concerning the model instances to be 
solved. 



TABLE 9-1 Problem Dimensions 





n 


m 


Pgub/% 


Ppn/% 


Pc;n/% 


NZEL 


GTE 


6,624 


960 


909/95 


909/95 


922/96 


58 


INVEST 


11,989 


1,338 


941/70 


1,101/82 


1,168/87 


33 


TANKER 


7,598 


83 


32/39 


32/39 


66/80 


31 


HFDF 


10,548 


61 


31/51 


31/51 


32/52 


189 


GAS PN A 


27,884 


6,848 


4,345/63 


5,934/87 


5,976/87 


37 


GAS PN C 


15,362 


3,794 


2,658/70 


3,418/90 


3,420/90 


20 


GAS PN E 


5,102 


1,184 


434/37 


877/74 


883/75 


7 


GAS GN A 


27,884 


6,848 


4,484/65 


5,142/75 


5,976/87 


37 


GAS GN C 


15,362 


3,794 


2,664/70 


3,084/81 


3,420/90 


20 


KELLOGG 2 


17,841 


3,818 


1,265/33 


2,578/68 


2,596/68 


35 


KELLOGG 3 


27,490 


5,727 


2,295/40 


3,867/68 


3,892/68 


54 


KELLOGG 4 


37,139 


7,636 


2,428/32 


5,156/68 


5,188/68 


74 


KELLOGG 5 


46,788 


9,545 


3,388/36 


6,445/68 


6,484/68 


93 


ODS 1 


11,568 


3,023 


528/17 


540/18 


558/18 


21 


ODS 3 


23,993 


594 


490/82 


490/82 


490/82 


68 


TAM 5 


10,531 


438 


102/23 


132/30 


162/37 


94 


TAM 8 


6,104 


420 


118/28 


154/37 


196/47 


49 


TAM 12 


17,793 


629 


177/28 


231 /37 


294/47 


165 


PHOENIX 10 


6,884 


1,618 


206/13 


220/14 


1,153/71 


14 


PHOENIX 30 


17,212 


4,305 


293/07 


303/07 


3,604/84 


48 


EA6B 


12,247 


2,978 


1,610/54 


2,921 /98 


2,921/98 


16 


DEC 


14,518 


2,171 


677/31 


677/31 


1,088/50 


24 


AMMO 4H 


83,497 


13,963 


6,874/49 


12,892/92 


12,892/92 


129 


BUSCH 4 


7,997 


1,248 


649/52 


1,140/91 


1,148/92 


15 


BAR 


49,032 


7,446 


2,712/36 


4,575/61 


5,134/69 


102 



For each problem, the total number of structural variables is n, structural con- 
straints m, GUB rows found by the identification heuristic pcuo, pure network 
rows pp;v, generalized network rows Pon-, an d thousands of nonzero technological 
coefficients NZEL. For example, GTE has 58 thousand nonzero technological co- 
efficients for 6, 624 variables; GTE can be viewed as having 960 explicit and no 
factored rows, or as 909/95% GUB-factored and 960 — 909 = 51 explicit rows, as 
909 PN-factored rows, or as 922 GN-factored and 38 explicit rows. 
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Table 9-2 displays solution times for the sample problems. These CPU times exclude initial 
problem input, factored row identification, and final output — on average about 0.2 second per 
problem. 



TABLE 9-2 Solution Seconds 







AMDAHL 5995-700 


486/33MHz 


PC 






X-System 




IBM 


XS 


CPLEX 




none 


GUB 


PN 


GN 


OSL 


GN 


1.2 


2.0 


GTE 


9 


8 


8 


8 


43 


41 


65 


58 


INVEST 


4 


5 


4 


4 


23 


24 


52 


49 


TANKER 


8 


10 


10 


8 


6 


48 


8 


9 


HFDF 


100 


101 


101 


111 


40 


462 


581 


542 


GAS PN A 


2,376 


778 


50 


57 


291 


321 


1,714 


1,221 


GAS PN C 


4 


4 


3 


3 


79 


26 


185 


245 


GAS PN E 


72 


33 


4 


6 


13 


33 


165 


50 


GAS GN A 


1,115 


542 


294 


72 


312 


385 


3,532 


2,262 


GAS GN C 


4 


4 


3 


3 


71 


26 


186 


236 


KELLOGG 2 


3 


2 


2 


2 


69 


5 


219 


194 


KELLOGG 3 


5 


5 


5 


5 


144 


9 


513 


440 


KELLOGG 4 


48 


36 


26 


24 


500 


115 


1,274 


1,111 


KELLOGG 5 


1,122 


1,459 


320 


248 


1,210 


926 


2,615 


2,243 


ODS 1 


6 


3 


2 


5 


125 


22 


1,042 


414 


ODS 3 


6 


6 


6 


6 


23 


38 


58 


56 


TAM 5 


55 


60 


45 


40 


44 


231 


151 


86 


TAM 8 


24 


14 


14 


17 


21 


69 


77 


71 


TAM 12 


108 


112 


88 


106 


101 


312 


416 


378 


PHOENIX 10 


4 


2 


2 


2 


20 


10 


84 


73 


PHOENIX 30 


41 


25 


26 


9 


335 


69 


1,171 


1,239 


EA6B 


75 


33 


9 


9 


45 


44 


177 


244 


DEC 


189 


32 


42 


80 


71 


93 


317 


293 


AMMO 4H 


42 


41 


35 


46 


1,000 


277 


1,762 


1,597 


BUSCH 4 


5 


5 


4 


4 


26 


34 


55 


46 


BAR 


290 


212 


106 


114 


382 


480 


1,366 


1,222 



An AMDAHL 5995-700 running under IBM VM/CMS/XA with IBM VS FOR- 
TRAN 2.3.0 is used to render performance in CPU-seconds accurate to the precision 
shown for the basic “XS”-system using no factorization, compared with dynamic 
factorizations of “XS(GUB)”, pure network “XS(PN)”, and generalized network 
“XS(GN)” rows. “IBM-OSL” shows primal simplex performance on the same com- 
puter of the IBM Optimization Subroutine Library, Release 2 [1991]. “486/33” 

shows the clock-time performance of “XS(GN)” on a microcomputer (33 MHz Intel 
486 with 32 MB RAM) followed by that of “CPLEX” (Version 1.2) [1990] and of 
“CP LEX” (Version 2.0) [1992] on the same microcomputer. 



The original formulations of most of the test problems are strongly influenced by the mod- 
eler’s solution strategy. For instance, TANKER is endowed with a GUB structure which places 
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every binary variable in an associated set from which only one member can be chosen; this 
maximal CUB set is also sought for its tendency to yield nearly- integer solutions to the linear 
program. AMMO 4H is a multicommodity capacitated transshipment problem and thus is best 
suited to a pure network factorization; it was originally solved by dual decomposition rendering 
pure network sub-problems. PHOENIX is a multicommodity equipment replacement model 
closely following the generalized network factorization paradigm. 

The row structures in Table 9-2 have been found without modeler help by our automatic 
identification heuristic, yet they corroborate the modelers 5 intentions. TANKER reveals the 
same pure network rows as the GUB rows, or a generalized network constructed by identifying 
one additional row to be paired with each GUB row. PHOENIX exhibits a dominant structure 
that is clearly a generalized network. 

One would anticipate the factorization exploiting the dominant row structure to win com- 
putation tests. This is wrong more often than right. Table 9-2 shows that the more general 
factorization dominates the less general, with few exceptions: GTE’s relatively large GUB set, 
and the large pure networks in GAS PN A, GAS PN E, and AMMO 4H seem to satisfy our 
prior bias toward model-dictated factorizations. 

Table 9-2 also suggests that our myopic use of heuristics to automatically identify factored 
structure has its pitfalls. In a number of problem instances, we identify significantly larger fac- 
tored sets with the more general factorizations, yet we enjoy little improvement in computation 
times (e.g., INVEST, ODS 3, TAM 8). This suggests that the “quality” of a row factorization 
is not completely specified by its size. 

We have pursued this notion by inviting some of the modelers to guide our identification 
heuristic to precisely the row sets they intended. Some of the results have been striking: Wood 
[1989] reports significant improvements for problems in the GAS system. 

A few of our corresponding modelers have had the opportunity to build models from scratch 
with a particular factorization in mind. This admits model coercion and a wide range of well- 
known reformulation methods which we think can materially change both the size and quality 
of the result. Their early reports show promise. Among the models discussed here, the GAS 
and KELLOGG systems have been subsequently reformulated to enhance generalized networks, 
GTE has been re-engineered to further accentuate its dominant GUB set, and TANKER-like 
and many ODS models have been moved to a micro-computer; all these models are now larger, 
but much easier to solve. 

It is surprising and encouraging that the transition to more general factorizations seldom 
degrades performance much, even when few additional factored rows are won by the increased 
generality. This contradicts popular folklore that the more general factorizations demand sub- 
stantial, if not overwhelming increases in the resulting sizes of the factored structures. In 
fact, computational testing reported by others has usually been limited to models in which 
the number of explicit rows is in the range of one to twenty (e.g., Chen and Saigal [1977], 
Glover, Karney, Klingman and Russell [1978], Glover and Klingman [1981]). Our results are 
all the more remarkable given the lack of guidance from the modeler for the “intended” row 
factorization. 
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Table 9-3 shows the maximum size of the A X \ in terms of its nonzero elements. 

TABLE 9-3 Maximum Number of Elements in Explicit Transformation Kernel 



GTE 
INVEST 
TANKER 
HFDF 
GAS PN A 
GAS PN C 
GAS PN E 
GAS GN A 
GAS GN C 
KELLOGG 2 
KELLOGG 3 
KELLOGG 4 
KELLOGG 5 
ODS 1 
ODS 3 
TAM 5 
TAM 8 
TAM 12 
PHOENIX 10 
PHOENIX 30 
EA6B 
DEC 

AMMO 4H 
BUSCH 4 
BAR 



XS XS(GUB) 



13 


0 


9 


1 


1 


1 


3 


1 


1,550 


891 


18 


5 


263 


110 


1,470 


758 


19 


7 


6 


2 


9 


4 


79 


41 


1,299 


1,015 


9 


1 


0 


0 


28 


22 


20 


15 


38 


42 


32 


21 


169 


157 


1,155 


270 


218 


39 


235 


92 


10 


5 


290 


163 



XS(PN) XS(GN) 



0 


0 


1 


1 


1 


0 


1 


1 


30 


54 


0 


0 


7 


5 


340 


59 


1 


0 


0 


0 


0 


0 


10 


11 


138 


128 


1 


5 


0 


0 


15 


18 


8 


10 


16 


21 


21 


1 


185 


1 


0 


0 


49 


46 


0 


0 


0 


0 


70 


41 



The number of nonzero elements (in nearest thousands) in the explicit transforma- 
tion kernel A ^ gives some indication of how much information is not captured by 
factorization, as well as an idea of relative storage requirements. 



We see that the maximum size of the explicit transformation kernel tends to decrease as the 
generality of the factorization increases. Recalling the definition of the explicit transformation 
kernel, 



An 1 = (E u -E l2 F 22 l F 2l y l , 

this trend is as we would expect. Each potentially binding explicit row which can be converted 
to a factored row reduces the likely size of . Also, the density of the term -EuF^ Fix 
generally increases with the density of ^22 ■ For Fxi fc-by-fc, the number of nonzeros in F 2 2 
for the GUB factorization is k , for pure networks perhaps as large y, and for generalized 
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networks as large as k 2 . There are some exceptions to this trend in Table 9-3, especially in 
model instances in which the size of the (GN) factored row set is not significantly larger than 
that of (PN). This is because for a given factored kernel F 22 , the exact (PN) representation of 
Foo is generally more sparse than that of the less exact floating-point representation of (GN). 
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It is usually the case that many constraints are not binding at optimality, as can be seen in 
Table 9-4. 



TABLE 9-4 Binding Explicit Constraints at Optimality 





Binding/% 


GUB/% 


PN/% 


GN/% 


GTE 


554/58 


20/02 


20/02 


18/02 


INVEST 


763/57 


201/15 


195/15 


162/12 


TANKER 


51/61 


60/72 


39/47 


9/11 


HFDF 


58/95 


29/48 


29/48 


28/46 


GAS PN A 


3,068/45 


1,953/29 


299/04 


349/05 


GAS PN C 


2,324/61 


850/22 


68/02 


90/02 


GAS PN E 


901/76 


573/48 


90/08 


79/07 


GAS GN A 


3,064/45 


1,854/27 


1,155/17 


359/05 


GAS GN C 


2,323/61 


861/23 


402/11 


89/02 


KELLOGG 2 


1,950/51 


1,015/27 


92/02 


118/03 


KELLOGG 3 


2,942/51 


1,368/24 


148/03 


183/03 


KELLOGG 4 


4,136/54 


2,491/33 


391 /05 


452/06 


KELLOGG 5 


5,270/55 


2,305/24 


592/06 


617/06 


ODS 1 


361/12 


83/03 


76/03 


117/04 


ODS 3 


410/69 


0/00 


0/00 


0/00 


TAM 5 


264/60 


227/52 


168/38 


186/42 


TAM 8 


279/66 


240/57 


142/34 


175/42 


TAM 12 


412/66 


352/53 


205/33 


251/40 


PHOENIX 10 


1,098/68 


1,096/68 


1,090/67 


80/05 


PHOENIX 30 


3,298/77 


3,236/75 


3,239/75 


110/03 


EA6B 


2,939/99 


1,334/45 


26/01 


27/01 


DEC 


1,328/61 


736/34 


726/33 


421/19 


AMMO 4H 


6,686/46 


3,153/22 


13/00 


14/00 


BUSCH 4 


840/67 


481 /39 


5/00 


8/01 


BAR 


4,687/63 


3,250/44 


1,895/25 


1,450/19 



Not all constraints are binding at optimality. The first column lists the number 
of binding explicit constraints and expresses this as a percentage of all constraints; 
the following columns display the number of binding explicit constraints and their 
corresponding percentages under the alternate factorizations. 



A distinguishing feature of dynamic factorization is the ability to limit attention to binding 
constraints, handling binding factored constraints with great efficiency, and working with a 
relatively small number of binding explicit constraints. Explicit binding constraints on the 
order of a few thousand, or less, are quite manageable. This is well beyond the size of previously 
reported implementations. 
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10 Conclusions 



Previous research by others generally suggests that specialized algorithms such as those pre- 
sented here are useful only when the factored structure completely dominates. There are even 
reports of algorithms for solving problems having a single unfactored (explicit) constraint (Hultz 
and Klingman [1978], Klingman and Russell [1978]). When implementations have been re- 
ported, problem suites have been limited to instances having a very small number of explicit 
constraints, typically in the range from one to twenty (Chen and Saigal [1977], Glover, Karney, 
Klingman and Russell [1978], Glover and Klingman [1981]). The consensus seems to be that 
such algorithms are quite delicate, and deserve to be viewed as specialized algorithms, useful 
only for solving very special problem instances. 

We refute this view. Dynamic factorization is competitive with commercial-quality opti- 
mization systems on every model instance we have tested. 

The development here stresses the similarities among the algorithms and the natural exten- 
sions leading from one to the next. This is in contrast to the development reported for similar, 
non-dynamic algorithms (e.g., Dantzig and Van Slyke [1967], Klingman and Russell [1978], 
and Hultz and Klingman [1978]) in which the specifics of the individual algorithm obscure the 
generality of the approach. The conceptual difference between our algorithms is seen to be 
largely isolated to the structure of a single algebraic entity, the factored kernel. By abstracting 
the structure of the factored kernel and concentrating on the general algorithm design, the 
versatility and flexibility of this approach is clarified. 

The algorithmic development leads directly to an implementation. The resulting software 
suite exhibits a “single system image”. The modularity of the algorithm allows the definition 
of an “abstract data type” (see, e.g., Aho, Hopcroft and Ullman [1974]) which isolates the data 
structures and update procedures for the factored kernel from the rest of the implementation. 
Each factorization is seamlessly integrated within the system design. 

The early 1980s produced a great deal of research on automatic identification of special 
structure in LP models (see, e.g,, Gunawardane and Schrage [1977], Glover [1980], Schrage 
[1981], Brown, McBride and Wood [1985], and Bixby and Fourer [1986]). We have incorporated 
the most useful of these ideas into our implementation, and we have what we believe to be the 
first complete implementation which supports automatic identification of factored row sets. This 
capability may be used to identify new factored structure or to validate or augment a modeler- 
provided recommendation. When faced with the choice of either solving an unfactored model 
instance or automatically identifying a factored structure and then using the corresponding 
solver, our results show that the latter is nearly always to be preferred. Modelers have conducted 
extensive additional computational experimentation with the X-System not reported in this 
paper. These results suggest that in addition to the quantity of factored rows, the quality of 
these rows influences the performance of factorization algorithms. While not well understood, 
it is clear that the myopic approach of our heuristics is no substitute for the modeler’s guidance 
in identifying factored structure. 

Processing networks (Koene [1982]) are network models which allow proportional flow re- 
strictions on the arcs entering or leaving some nodes. One formulation of such a model results 
in a pure or generalized network structure with a set of complicating columns. Chen and En- 
gquist [1986] propose a primal partitioning algorithm for solving processing network problems. 
An alternate formulation yields a pure or generalized network structure with complicating rows: 
this is precisely the structure dealt with here. 

The multicommodity capacitated transshipment problem (MCTP) has been the subject of 
much research over the years, and a number of specialized algorithms have been proposed to 
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solve it (see, e.g., Assad [1978] or Kennington [1978]). The usual MCTP formulation is a pure 
network which each commodity uses independently in its own flow model, but with side con- 
straints on the total common flow of all commodities over some of the network arcs. The side 
constraints form a GUB row set, while the rest of MCTP forms a pure network; either view 
might be preferred depending upon size of the common network, the number of side constraints, 
and the number of commodities. In our experience, the network factorization usually domi- 
nates the GUB factorization, and the pure network factorization presented here is a powerful 
technique for solving MCTPs. As an experiment, we customized our (PN) implementation for 
MCTP to exploit the special structure of the explicit side constraints. This highly-specialized 
implementation performed no better on AMMO 4H, and we now believe that this would be 
true for most MCTPs. 

There are problems which would exhibit a large factored row structure if not for a set of 
complicating columns (e.g., see Brown, McBride and Wood [1985]). One would expect the 
structure of the factored kernel to be dominated by that of the predominant row structure, 
with only occasional complications due to the exceptional columns. One might allow for this 
exceptional structure in the factored kernel by identifying it “on-the-fly” as the algorithm 
progresses, and preserving the sanctity of the core factorization. Though conceptually simple, 
some iterations of this algorithm would border on the spectacular. This approach may be 
thought of as a hybrid between the dynamic factorization developed here and dynamic basis 
triangulation methods (see, e.g., Hellerman and Rarick [1971] and [1972], Saunders [1976] and 
McBride [1980]). 

Dynamic extrinsic factorization is subsumed by the algorithms presented in this paper if 
we activate functions in the update analogous to the secondary exchanges now employed. Es- 
sentially all that has to be done is ensure that successive factored components retain their 
stipulated special structure. We speculate that this will work best in cases where model struc- 
ture is amenable, and quite likely will require some model-specific customization to perform 
well on difficult models. We have limited our experimentation to those static extrinsic cases 
which we believe to be most generally useful. 
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